home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 09 - 1993 / 09.10 Oct 93 / Simpson's Rule / Auxiliary.p next >
Encoding:
Text File  |  1993-07-25  |  3.0 KB  |  88 lines  |  [TEXT/PJMM]
  1. { © 1993  Marek Hajek. Hajek 's Solutions. Created 6-12-93 }
  2. {--This example illustrates Integral Computations }
  3.  
  4. UNIT Auxiliary;
  5. INTERFACE
  6.     USES
  7.         Sane;
  8.  
  9.     FUNCTION ComputeIntegral (lowerLimit, upperLimit: Extended;
  10.                                     partitionCount: LongInt;
  11.                                     FUNCTION IntegrandFunction (partitionCoordinate: Extended): Extended): Extended;
  12.  
  13. (* You have to define the Integrand function *)
  14.     FUNCTION IntegrandFunction (partitionCoordinate: Extended): Extended;
  15. (* Computes the head of the equation *)
  16.     FUNCTION ComputeHead (lowerLimit, upperLimit: Extended;
  17.                                     partitionCount: LongInt): Extended;
  18. (* Computes first and last of the equation *)
  19.     FUNCTION FirstAndLast (lowerLimit, upperLimit: Extended;
  20.                                     FUNCTION IntegrandFunction (partitionCoordinate: Extended): Extended): Extended;
  21.  
  22. IMPLEMENTATION
  23.  
  24. {--------------------------IntegrandFunction----------------------------}
  25.     FUNCTION IntegrandFunction (partitionCoordinate: Extended): Extended;
  26. (* This function is the one inside the integral *)
  27. (*Translate it into pascal and define it here *)
  28.  
  29. { This functions computes -> √(X * X * X * X +1 }
  30.     BEGIN
  31.         IntegrandFunction := SQRT(XpwrI(partitionCoordinate, 4) + 1);
  32.     END;
  33.  
  34. {--------------------------ComputeHead----------------------------}
  35.     FUNCTION ComputeHead (lowerLimit, upperLimit: Extended;
  36.                                     partitionCount: LongInt): Extended;
  37. (*  Computes the first part of the integral equation, the Head  *)
  38. (*  Corresponds to (b - a)/(3*n)  *)
  39.  
  40.     BEGIN
  41.         ComputeHead := (upperLimit - lowerLimit) / (3 * partitionCount);
  42.     END;
  43.  
  44. {--------------------------FirstAndLast----------------------------}
  45.     FUNCTION FirstAndLast (lowerLimit, upperLimit: Extended;
  46.                                     FUNCTION IntegrandFunction (partitionCoordinate: Extended): Extended): Extended;
  47. (*  Computes the second part of the integral, the FIRST/LAST  *)
  48. (*  Corresponds to [f(X0) + f(Xn)  *)
  49.  
  50.     BEGIN
  51.         FirstAndLast := IntegrandFunction(lowerLimit) + IntegrandFunction(upperLimit);
  52.     END;
  53.  
  54. {--------------------------ComputeIntegral----------------------------}
  55.     FUNCTION ComputeIntegral (lowerLimit, upperLimit: Extended;
  56.                                     partitionCount: LongInt;
  57.                                     FUNCTION IntegrandFunction (partitionCoordinate: Extended): Extended): Extended;
  58.         VAR
  59.             result: Extended;
  60.             head: Extended;
  61.             partitionIncrement: Extended;
  62.             partitionCoordinate: Extended;
  63.             index: LongInt;
  64.  
  65.     BEGIN
  66.         head := ComputeHead(lowerLimit, upperLimit, partitionCount);
  67.         result := FirstAndLast(lowerLimit, upperLimit, IntegrandFunction);
  68.  
  69.         partitionIncrement := (upperLimit - lowerLimit) / partitionCount;
  70.         partitionCoordinate := lowerLimit;
  71.  
  72.         FOR index := 1 TO partitionCount - 1 DO
  73.             BEGIN
  74. (* Partition coordinate coresponds to X0, X1, X2,.....Xn *)
  75.                 partitionCoordinate := partitionCoordinate + partitionIncrement;
  76.  
  77. (* Odd index means compute 4 * f(x), even index means compute 2 * f(x)  *)
  78.                 IF Odd(index) THEN
  79.                     result := result + 4 * IntegrandFunction(partitionCoordinate)
  80.                 ELSE
  81.                     result := result + 2 * IntegrandFunction(partitionCoordinate)
  82.  
  83.             END;  (* FOR ... *)
  84.  
  85.         ComputeIntegral := head * result;
  86.     END;
  87.  
  88. END.